Loading libraries and set directories

library(dplyr)
library(apeglm)
library(org.Hs.eg.db)
library(Seurat)
library(sctransform)
library(SeuratObject)
library(DESeq2)
library(enrichplot)
library(clusterProfiler)
library(AnnotationDbi)
library(SingleCellExperiment)
library(muscat)
library(ggplot2)
library(gprofiler2)
library(ggrepel)
library(msigdbr)


MEDIApath="/home/mclaudia/MEDIA Project/"
MEDIAgraph="/home/mclaudia/MEDIA Project/Codici_finali/Paper_figures/"
MEDIAsave="/home/mclaudia/MEDIA Project/Codici_finali/"
wd="/home/mclaudia/MEDIA Project/GSE152805_RAW_scell/"

SingleCell Dataset from Chou et al., 2020 (doi:10.1038/s41598-020-67730-y): Dataset 1

Build merged Seurat Object of the three patients: s113, s116, s118

MT: medial region of the tibia (affected cartilage)

OLT: outer lateral region of the tibia (healthy cartilage)

MT1.data <-Read10X(data.dir =paste0(wd,"MT_113"))
MT1<- CreateSeuratObject(counts = MT1.data, 
                         project = "MT1", min.cells = 300, min.features = 500)

MT2.data <-Read10X(data.dir = paste0(wd,"MT_116"))
MT2<- CreateSeuratObject(counts = MT2.data, 
                         project = "MT2", min.cells = 300, min.features = 500)

MT3.data <-Read10X(data.dir = paste0(wd,"MT_118"))
MT3<- CreateSeuratObject(counts = MT3.data, 
                         project = "MT3", min.cells = 300, min.features = 500)

olt1.data <-Read10X(data.dir = paste0(wd,"olt_113"))
olt1<- CreateSeuratObject(counts = olt1.data, 
                          project = "olt1", min.cells = 300, min.features = 500)

olt2.data <-Read10X(data.dir = paste0(wd,"olt_116"))
olt2<- CreateSeuratObject(counts = olt2.data, 
                          project = "olt2", min.cells = 300, min.features = 500)

olt3.data <-Read10X(data.dir = paste0(wd,"olt_118"))
olt3<- CreateSeuratObject(counts = olt3.data, 
                          project = "olt3", min.cells = 300, min.features = 500)

tot113 <- merge(MT1, y = olt1, add.cell.ids = c("MT_s113", "Olt_s113"), project = "tot113")
length(Cells(tot113))

tot116 <- merge(MT2, y = olt2, add.cell.ids = c("MT_s116", "Olt_s116"), project = "tot116")
length(Cells(tot116))

tot118 <- merge(MT3, y = olt3, add.cell.ids = c("MT_s118", "Olt_s118"), project = "tot118")
length(Cells(tot118))

pt_list <-list(tot113,tot116,tot118)
names(pt_list) <-c("tot113","tot116","tot118")

Pre-processing

for( i in 1:length(pt_list)){
  
  counts <- GetAssayData(object = pt_list[[i]], slot = "counts")
  
  nonzero <- counts > 0
  keep_genes <- Matrix::rowSums(nonzero) >= 100
  filtered_counts <- counts[keep_genes, ]
  filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = pt_list[[i]]@meta.data)
  
  tmp <-filtered_seurat
  
  tmp[["percent.mt"]] <- PercentageFeatureSet(tmp , pattern = "^MT-")

  VlnPlot(tmp, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  
  plot1 <- FeatureScatter(tmp, feature1 = "nCount_RNA", feature2 = "percent.mt")
  plot2 <- FeatureScatter(tmp, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
 # plot1 + plot2, visualize to set nFeature_RNA, nCount_RNA and percent.mt thresholds
  
  if(i == 1) tmp <- subset(tmp, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 4 & nCount_RNA < 50000)
  if(i == 2) tmp <- subset(tmp, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 5 & nCount_RNA < 50000)
  if(i == 3) tmp <- subset(tmp, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 3 & nCount_RNA < 40000)
  
  tmp<- NormalizeData(tmp)
  tmp1<- FindVariableFeatures(tmp, selection.method = "vst", nfeatures = 8500)
  pt_list[[i]] <-tmp1
}

Integrate patient data

features <- SelectIntegrationFeatures(pt_list, nfeatures = 8500) 

anchors <- FindIntegrationAnchors(pt_list, anchor.features=features)

OAtot <-IntegrateData(anchors, new.assay.name = "integrate")

Examine and visualize PCA

all.genes <-  rownames(OAtot)

OAtot <- ScaleData(OAtot, features = all.genes)
OAtot <- RunPCA(OAtot, features = VariableFeatures(object = OAtot))
print(OAtot[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1 
## Positive:  S100A4, COL1A2, CRTAC1, HTRA1, COL3A1 
## Negative:  FRZB, CFH, JUNB, ACAN, COL9A3 
## PC_ 2 
## Positive:  CHI3L2, NNMT, CHI3L1, CLU, RPL13 
## Negative:  FGFBP2, S100A1, SCRG1, C2orf82, CHAD 
## PC_ 3 
## Positive:  COMP, OMD, PRELP, CILP, SMOC2 
## Negative:  LMNA, COL1A1, VIM, NBL1, B2M 
## PC_ 4 
## Positive:  IBSP, COL10A1, SPP1, MT-ND3, SOD3 
## Negative:  COL2A1, EEF1A1, RPL3, CILP2, RPS2 
## PC_ 5 
## Positive:  RPL23A, RPS16, RPL37, RPL37A, RPLP2 
## Negative:  CLU, ITM2B, LUM, DCN, PDIA3

Run UMAP: 10 PC chosen

VizDimLoadings(OAtot, dims = 1:2, reduction = "pca")

DimPlot(OAtot, reduction = "pca")

DimHeatmap(OAtot, dims = 1:15, cells = 500, balanced = TRUE)

ElbowPlot(OAtot) ###9-12 PC

OAtot <- FindNeighbors(OAtot, dims = 1:10)
OAtot <- FindClusters(OAtot, resolution = 0.5) 
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 23752
## Number of edges: 690120
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8396
## Number of communities: 9
## Elapsed time: 12 seconds
OAtot <- RunUMAP(OAtot, dims = 1:10)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

Save UMAP coordinates for reproducibility

umapCoord <- Embeddings(object = OAtot[["umap"]])
save(umapCoord,file=paste0(MEDIAsave,"CoordUMAP.RData"))

Visualize UMAP for: clusterization, class, single patient

load(paste0(MEDIAsave,"CoordUMAP.RData"))
OAtot@meta.data$type <-sapply(strsplit(rownames(OAtot@meta.data),"_"), function(x) x[1])
OAtot@meta.data$patient <-sapply(strsplit(rownames(OAtot@meta.data),"_"), function(x) x[2])

OAtot@reductions[["umap"]]@cell.embeddings  <- umapCoord

DimPlot(OAtot, reduction = "umap", label = TRUE, label.size = 7)+
   ggtitle("Clustering")+
   theme(plot.title = element_text(hjust = 0.5))

DimPlot(OAtot, reduction = "umap", group.by = "type", label = TRUE, 
        label.size = 7,cols= c("turquoise","orangered"))+
        ggtitle("Class")

DimPlot(OAtot, reduction = "umap", group.by = "patient")+
        ggtitle("Patient")

Pseudo-bulk singleCell data trasformation

ei <-OAtot@meta.data[,c(1,6,7,8)]

head(ei) 
##                            orig.ident seurat_clusters type patient
## MT_s113_AAACCTGGTAGCACGA-1        MT1               2   MT    s113
## MT_s113_AAACGGGAGGCAGGTT-1        MT1               2   MT    s113
## MT_s113_AAACGGGCAAATACAG-1        MT1               2   MT    s113
## MT_s113_AAACGGGCAAGAGTCG-1        MT1               0   MT    s113
## MT_s113_AAACGGGCACCGCTAG-1        MT1               0   MT    s113
## MT_s113_AAACGGGCAGTCTTCC-1        MT1               0   MT    s113
ei$clus <-0  #merge in unique

OA_sce <-as.SingleCellExperiment(OAtot, assay = "RNA") 
OA_sce$clus <-0

(OA_sce2 <- prepSCE(OA_sce, 
                    kid = "seurat_clusters",
                    gid = "type",  
                    sid = "orig.ident",   
                    drop = FALSE))  
## class: SingleCellExperiment 
## dim: 10060 23752 
## metadata(1): experiment_info
## assays(2): counts logcounts
## rownames(10060): AP006222.2 SAMD11 ... LCA5L LINC00205
## rowData names(0):
## colnames(23752): MT_s113_AAACCTGGTAGCACGA-1 MT_s113_AAACGGGAGGCAGGTT-1
##   ... Olt_s118_TTTGTCAGTGCACGAA-1 Olt_s118_TTTGTCATCCTAAGTG-1
## colData names(10): cluster_id sample_id ... ident clus
## reducedDimNames(2): PCA UMAP
## mainExpName: RNA
## altExpNames(0):
OA_sce2$cond <-paste0(OA_sce2$sample_id,"_",OA_sce2$patient)

experiment_info <- data.frame(sample_id=as.factor(c("MT1","MT2","MT3","olt1","olt2","olt3")),
                              cond=as.factor(c("MT1_s113 ","MT2_s116",
                                               "MT3_s118","olt1_s113","olt2_s116","olt3_s118")),
                              patient=as.factor(rep(c("s113","s116","s118"),2)),
                              group_id=as.factor(c("MT","MT","MT","Olt","Olt","Olt")),
                              n_cells=as.numeric(table(ei$orig.ident)))


OA_sce2@metadata <- list(experiment_info)
names(OA_sce2@metadata) <-"experiment_info"


table(OA_sce2$cluster_id,OA_sce2$sample_id)
##    
##      MT1  MT2  MT3 olt1 olt2 olt3
##   0 1733  903 2477  223  321  325
##   1   64   27   62 2260 1135 2037
##   2  897  442 1189  101   67   70
##   3    2   13   81 1350  933  378
##   4    0    0   12  858  895  633
##   5   18  134   42 1115  259  474
##   6  243  693  377   16    5   66
##   7   51  486    4    2    0   20
##   8   22    0   22  155   46   14
pb <- aggregateData(OA_sce2, assay = "counts", fun = "sum")
colData(pb)
## DataFrame with 6 rows and 4 columns
##      group_id     patient      clus        cond
##      <factor> <character> <numeric> <character>
## MT1       MT         s113         0    MT1_s113
## MT2       MT         s116         0    MT2_s116
## MT3       MT         s118         0    MT3_s118
## olt1      Olt        s113         0   olt1_s113
## olt2      Olt        s116         0   olt2_s116
## olt3      Olt        s118         0   olt3_s118
plotPCA(OA_sce2, colour_by = "group_id")

pbMDS(pb)

colData(pb)$patient <-as.factor(colData(pb)$patient)
colData(pb)$cond <-as.factor(colData(pb)$cond)

ei <- metadata(pb)$experiment_info

pb$group_id <-relevel(pb$group_id,ref="Olt")
ei$group_id <-relevel(ei$group_id,ref="Olt")


ass <-matrix(0, nrow = nrow(pb@int_elementMetadata), ncol = 6)
colnames(ass) <-c("MT1","MT2","MT3","olt1","olt2","olt3")

for(i in c(1:9)){
  as <-as.data.frame(as.matrix(pb@assays@data@listData[[i]]))
  ass <-ass+as
}

ass <- ass+1  #add dummy value = 1

cd <- colData(pb)
y <- DESeqDataSetFromMatrix(ass, cd, design=~patient+group_id)

dds <- DESeq(y)
resultsNames(dds)
## [1] "Intercept"            "patient_s116_vs_s113" "patient_s118_vs_s113"
## [4] "group_id_MT_vs_Olt"
res <-results(dds,alpha=.9,name="group_id_MT_vs_Olt",pAdjustMethod = "fdr")

df.res_sc <-lfcShrink(dds, coef="group_id_MT_vs_Olt",res =res,type = "apeglm")
df.res_sc <- as.data.frame(df.res_sc)
df.res_sc <-df.res_sc[order(df.res_sc$padj, decreasing = FALSE),]

df.res_sc[1:10,]
##         baseMean log2FoldChange      lfcSE       pvalue         padj
## COL1A1 12787.080      14.356062 0.83122461 1.311834e-77 1.319705e-73
## ABI3BP  9863.181       1.605177 0.09261043 2.269700e-68 1.141659e-64
## TGFBI  12216.317       2.863385 0.16669522 7.913243e-68 2.653574e-64
## TSPAN2 13968.174       1.685722 0.10144715 3.349151e-63 8.423114e-60
## PDLIM7  2412.373       2.232388 0.14327520 4.262484e-56 8.576118e-53
## CRLF1   3719.325      12.742026 0.81759672 2.333293e-52 3.912155e-49
## MT1G   77050.541       1.863951 0.12565511 5.386437e-51 7.741079e-48
## CRTAC1 46489.316       2.173891 0.14739598 1.417025e-50 1.781909e-47
## KLF10   6336.885      -1.184766 0.08453508 9.633670e-46 1.076830e-42
## CLIC3   1633.566      11.641215 0.81663338 4.616620e-44 4.644320e-41

Save files

save(df.res_sc,file=paste0(MEDIAsave,"DE_DESeq2_SCell.RData"))

Volcano Plot of DE genes

annot <-df.res_sc[,c(2,5)]
annot$gene_names <-rownames(df.res_sc)

annot$col <-"gray50"
annot[annot$log2FoldChange <= -9 & annot$padj <=0.05, "col"] <- "springgreen3"
annot[annot$log2FoldChange  >= 9 & annot$padj <=0.05, "col"] <- "red"

annot$label <- NA
annot[annot$log2FoldChange  >= 10 & annot$padj <=0.05, "label"] <- annot[annot$log2FoldChange  >= 10 & annot$padj <=0.05, "gene_names"]
annot[annot$log2FoldChange  <= -10 & annot$padj <=0.05, "label"] <- annot[annot$log2FoldChange   <= -10 & annot$padj <=0.05, "gene_names"]

ggplot(data=annot, aes(x=log2FoldChange, y=-log10(padj), color=col,label=label)) +
  geom_point(size=1) +
  geom_label_repel(size=4,
                   min.segment.length = 0,direction = "both",
                   max.overlaps =40,
                   color="black",
                   segment.linetype=2) +
  scale_shape_manual(values=c(8, 19))+
  scale_color_manual(values = c("red","gray50","springgreen3"),
                     breaks = c("red","gray50","springgreen3"),
                     labels = c("UP","notDE","DOWN")) +
  labs(color="Legend")+
  xlab("log2(FC)")+
  ylab("-log10(FDR)")+
  theme_bw()+
  ggtitle("Volcano plot of DE genes - Dataset 1")+
  geom_vline(xintercept=c(-9, 9), col="gray70", linetype=2)+
  geom_hline(yintercept=-log10(0.05), col="gray70", linetype=2)

Enrichment with GSEA: BP and REACTOME

hs=get("org.Hs.eg.db")

genes1 <- annot[, c("log2FoldChange","gene_names")]
genes <- genes1$log2FoldChange
names(genes) <-genes1$gene_names

genes <-sort(genes, decreasing = TRUE)

gsebp <- gseGO(geneList=genes, 
               ont ="BP", 
               keyType = "SYMBOL", 
               pvalueCutoff = 0.01,
               eps = 0,
               verbose = TRUE, 
               OrgDb = hs, 
               pAdjustMethod = "fdr")

dotplot(gsebp, showCategory = 20, x="NES",split=".sign") + facet_grid(.~.sign)

rc <-msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:REACTOME")
sig <-rc[,c(3,4)]
sig$gs_name <-gsub("^REACTOME_","",sig$gs_name)

gsreac <- clusterProfiler::GSEA(genes,TERM2GENE = sig,
                             eps = 0,
                             verbose = TRUE, 
                             minGSSize =5, pAdjustMethod = "fdr",
                             pvalueCutoff =0.1)
dotplot(gsreac, showCategory = 10, x="NES",split=".sign") + facet_grid(.~.sign)+
  theme(axis.text.y = element_text(size = 8))

Session Info

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] msigdbr_7.5.1               ggrepel_0.9.3              
##  [3] gprofiler2_0.2.2            ggplot2_3.4.3              
##  [5] muscat_1.10.1               SingleCellExperiment_1.18.1
##  [7] clusterProfiler_4.4.4       enrichplot_1.21.0          
##  [9] DESeq2_1.36.0               SummarizedExperiment_1.26.1
## [11] MatrixGenerics_1.8.1        matrixStats_1.0.0          
## [13] GenomicRanges_1.48.0        GenomeInfoDb_1.32.4        
## [15] sctransform_0.3.5           SeuratObject_4.1.3         
## [17] Seurat_4.3.0.1              org.Hs.eg.db_3.15.0        
## [19] AnnotationDbi_1.58.0        IRanges_2.30.1             
## [21] S4Vectors_0.34.0            Biobase_2.56.0             
## [23] BiocGenerics_0.42.0         apeglm_1.18.0              
## [25] dplyr_1.1.3                
## 
## loaded via a namespace (and not attached):
##   [1] rsvd_1.0.5                ica_1.0-3                
##   [3] foreach_1.5.2             lmtest_0.9-40            
##   [5] crayon_1.5.2              rbibutils_2.2.13         
##   [7] MASS_7.3-58.3             nlme_3.1-162             
##   [9] backports_1.4.1           GOSemSim_2.22.0          
##  [11] rlang_1.1.1               XVector_0.36.0           
##  [13] ROCR_1.0-11               irlba_2.3.5.1            
##  [15] nloptr_2.0.3              limma_3.52.4             
##  [17] scater_1.24.0             BiocParallel_1.30.4      
##  [19] rjson_0.2.21              bit64_4.0.5              
##  [21] glue_1.6.2                pbkrtest_0.5.2           
##  [23] parallel_4.2.1            vipor_0.4.5              
##  [25] spatstat.sparse_3.0-2     DOSE_3.22.1              
##  [27] spatstat.geom_3.2-1       tidyselect_1.2.0         
##  [29] fitdistrplus_1.1-11       variancePartition_1.28.9 
##  [31] XML_3.99-0.14             tidyr_1.3.0              
##  [33] zoo_1.8-12                xtable_1.8-4             
##  [35] magrittr_2.0.3            evaluate_0.21            
##  [37] Rdpack_2.4                scuttle_1.8.4            
##  [39] cli_3.6.1                 zlibbioc_1.42.0          
##  [41] rstudioapi_0.15.0         miniUI_0.1.1.1           
##  [43] sp_2.0-0                  bslib_0.5.1              
##  [45] fastmatch_1.1-4           aod_1.3.2                
##  [47] EnvStats_2.7.0            treeio_1.20.2            
##  [49] shiny_1.7.5               BiocSingular_1.14.0      
##  [51] xfun_0.39                 clue_0.3-64              
##  [53] cluster_2.1.4             caTools_1.18.2           
##  [55] tidygraph_1.2.3           KEGGREST_1.36.3          
##  [57] tibble_3.2.1              ape_5.7-1                
##  [59] listenv_0.9.0             Biostrings_2.64.1        
##  [61] png_0.1-8                 future_1.33.0            
##  [63] withr_2.5.0               bitops_1.0-7             
##  [65] ggforce_0.4.1             plyr_1.8.8               
##  [67] coda_0.19-4               pillar_1.9.0             
##  [69] gplots_3.1.3              GlobalOptions_0.1.2      
##  [71] cachem_1.0.8              multcomp_1.4-25          
##  [73] GetoptLong_1.0.5          DelayedMatrixStats_1.20.0
##  [75] vctrs_0.6.3               ellipsis_0.3.2           
##  [77] generics_0.1.3            tools_4.2.1              
##  [79] remaCor_0.0.16            beeswarm_0.4.0           
##  [81] munsell_0.5.0             tweenr_2.0.2             
##  [83] fgsea_1.22.0              emmeans_1.8.8            
##  [85] DelayedArray_0.22.0       fastmap_1.1.1            
##  [87] compiler_4.2.1            abind_1.4-5              
##  [89] httpuv_1.6.11             plotly_4.10.2            
##  [91] GenomeInfoDbData_1.2.8    gridExtra_2.3            
##  [93] glmmTMB_1.1.7             edgeR_3.38.4             
##  [95] lattice_0.21-8            deldir_1.0-9             
##  [97] utf8_1.2.3                later_1.3.1              
##  [99] jsonlite_1.8.7            scales_1.2.1             
## [101] ScaledMatrix_1.6.0        tidytree_0.4.2           
## [103] pbapply_1.7-2             sparseMatrixStats_1.10.0 
## [105] estimability_1.4.1        genefilter_1.78.0        
## [107] lazyeval_0.2.2            promises_1.2.1           
## [109] doParallel_1.0.17         goftest_1.2-3            
## [111] spatstat.utils_3.0-3      reticulate_1.31          
## [113] rmarkdown_2.24            sandwich_3.0-2           
## [115] cowplot_1.1.3             blme_1.0-5               
## [117] Rtsne_0.16                downloader_0.4           
## [119] uwot_0.1.16               igraph_1.5.0             
## [121] survival_3.5-5            numDeriv_2016.8-1.1      
## [123] yaml_2.3.7                htmltools_0.5.6          
## [125] memoise_2.0.1             locfit_1.5-9.8           
## [127] graphlayouts_1.0.0        viridisLite_0.4.2        
## [129] digest_0.6.32             RhpcBLASctl_0.23-42      
## [131] mime_0.12                 emdbook_1.3.13           
## [133] RSQLite_2.3.1             yulab.utils_0.0.6        
## [135] future.apply_1.11.0       data.table_1.14.8        
## [137] blob_1.2.4                splines_4.2.1            
## [139] labeling_0.4.3            RCurl_1.98-1.12          
## [141] broom_1.0.5               hms_1.1.3                
## [143] colorspace_2.1-0          ggbeeswarm_0.7.2         
## [145] shape_1.4.6               aplot_0.1.10             
## [147] sass_0.4.7                Rcpp_1.0.11              
## [149] RANN_2.6.1                mvtnorm_1.2-3            
## [151] circlize_0.4.15           fansi_1.0.4              
## [153] parallelly_1.36.0         R6_2.5.1                 
## [155] grid_4.2.1                ggridges_0.5.4           
## [157] lifecycle_1.0.3           minqa_1.2.5              
## [159] leiden_0.4.3              jquerylib_0.1.4          
## [161] DO.db_2.9                 Matrix_1.5-4.1           
## [163] qvalue_2.28.0             RcppAnnoy_0.0.21         
## [165] TH.data_1.1-2             RColorBrewer_1.1-3       
## [167] iterators_1.0.14          spatstat.explore_3.2-1   
## [169] TMB_1.9.4                 stringr_1.5.0            
## [171] htmlwidgets_1.6.2         beachmat_2.14.2          
## [173] polyclip_1.10-4           purrr_1.0.1              
## [175] shadowtext_0.1.2          gridGraphics_0.5-1       
## [177] ComplexHeatmap_2.12.1     globals_0.16.2           
## [179] patchwork_1.1.3           spatstat.random_3.1-5    
## [181] bdsmatrix_1.3-6           progressr_0.14.0         
## [183] codetools_0.2-19          GO.db_3.15.0             
## [185] gtools_3.9.4              prettyunits_1.1.1        
## [187] gtable_0.3.4              DBI_1.1.3                
## [189] ggfun_0.1.0               tensor_1.5               
## [191] httr_1.4.7                highr_0.10               
## [193] KernSmooth_2.23-20        stringi_1.7.12           
## [195] progress_1.2.2            reshape2_1.4.4           
## [197] farver_2.1.1              annotate_1.74.0          
## [199] viridis_0.6.3             ggtree_3.6.2             
## [201] bbmle_1.0.25              boot_1.3-28.1            
## [203] BiocNeighbors_1.16.0      lme4_1.1-34              
## [205] geneplotter_1.74.0        ggplotify_0.1.0          
## [207] scattermore_1.2           bit_4.0.5                
## [209] scatterpie_0.2.1          spatstat.data_3.0-1      
## [211] ggraph_2.1.0              pkgconfig_2.0.3          
## [213] babelgene_22.9            lmerTest_3.1-3           
## [215] knitr_1.43